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Various widely-used mean-field type theories for a dilute Bose gas are critically examined in 
the light of the recent discovery of Bose-Einstein condensation of atomic gases in a confined 
geometry. By numerically solving the mean-field equations within the framework of the Bogoli- 
ubov approximation both stationary non-uniform case and the vortex case under rotation in 
a cylindrically symmetric vessel are investigated. We obtain spatial structures of condensate, 
non-condensate, anomalous correlation. The low lying excitation spectra, the local density of 
states and the circulating current density in a vortex corresponding to various levels of mean-field 
theories are predicted. 
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§1. Introduction 

Much attention has been focused on many-body di- 
lute Bose systems since the discovery of Bose-Einstein 
condensation (BEC) in alkali atom gasjes-.ajt ultra-low 
temperatures arctund 0{fiK) ^ 0{nK) a na' The cur- 
rent experimentso' performed on atomic gases such as 
sodium and rubidium in confined geometries trapped by 
a magnetic field and/or optical method. Furious exper- 
imental and theoretical investigations are now ongoing 
and rapidly developing. Since these atom gases are di- 
lute Bose systems, allowing us to model these as weakly 
interacting many-body Bose systems, it is regarded as re- 
alization of ideal Bose-Einstein condensation. Although 
superfluid ^He has been considered as the first BEC, the 
mutual repulsive interaction is rather strong where the 
condensate fraction only ^--^10% of the total. There is 
little chance to directly apply microscopic theory to it. 

The microscopic theoretical work on BEC started with 
Bogoliubov J3' followed by several important progresses, 
such as Gross,© Pitaevskii,EP Iordanskii,EP and FetterB 
Thoz, treat an infinite system. These mean-field theo- 



rieal-l'' are particularly suited for treating an spatially 
non-uniform system which is the case of the present gaseous 
BEC systems trapped optically or magnetically in a re- 
stricted geometry. Therefore, the current theoretical stud- 
ieglil' mainly focus on examining these mean field theo- 
ries at various stages of the employed approximation for 
dilute Bose gases trapped by the usually harmonic po- 
tential in order to extract the spatial structures of the 
condensate and non-condensate and low-lying collective 
modes. Agreement between these mean field theories and 
experiments is generally fairly good so far, encouraging 
us to go along this line. 

The purposes of this paper are two-fold: One is to 
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compare various approximate treatments: (1) the Gross- 
Pitaevskii approximation (CP), which considers only the 
condensate, (2) the Bogoliubov approximation (BA), which 
takes into account the non-condensate too, but not self- 
consistent, (3) the Popov approximation (PA), which ne- 
glects the anomalous correlation and (4) the Hartree- 
Fock Bogoliubov approximation (HFB). (3) and (4) take 
into account the condensate and non-condensate self- 
consistently. We obtain a complete self-consistent solu- 
tion numerically for these four cases on an equal footing. 
This has not been done before. This is particularly true 
for HFB whose detailed study has not been performed. 
These various approximate theories reveal both merits 
and demerits. None is complete and satisfactory. Never- 
theless we may gain some insights to general properties 
of trapped BEC systems. 

The second purpose, which is the main purpose, is to 
know various properties of the quantized vortex struc- 
ture in rotating BEC systems since microscopic treat- 
ment of the vortex has beeii lacking so far except for a 
few attempts based on GPEj^ for present BEC systems. 
We are particularly interested in understanding (A) the 
spatial profiles of the condensate, non-condensate and 
anomalous correlation, (B) the dispersion relation of the 
excitation spectrum of the system, and (C) the low-lying 
excitations localized around a yortex core in detail, which 
is known as the Kelvin modeJld) These features should 
be compared with those in vortices in a supeixonduc- 
tor where extensive studies have been compiled.E-S This 
is particularly true for the localized, excitations near a 
vortex core studied by Caroli et alS3 The HFB theory 
just corresponds to the so-called Bogoliubov de Gcnnes 
theory (BdG), which is widely used to treat various spa- 
tially non-uniform superconductors havirig^ interface or 
vortex. In particular, Gygi and Schliitertj' have done 
a detailed investigation on an isolated vortex, providing 
a fairly complete and clear picture of microscopic struc- 
tures for quantized vortex in a s-wave superconductor. 
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We have been motivated by their impressive work and 
considered that a similar work should be done in a BEC 
system because vortex in BEC might be eventually ob- 
served and the vortex core structure and associated low- 
lying excitations may be probed experimentally. The 
characteristic core radius is estimated as the coherence 
length ~ 0.1/im which contrasted with a few A in super- 
fluid ^He. Thus there is a good chance to investigate the 
detailed core structure of the present BEC systems. In 
a superconductor scanning electron microscopy by Hess 
et al.t3 directly gives vivid spatially resolved images of 
the low-lying excitations near a vortex core. 

It turns out that although the full self-consistent HFB 
solution obtained here shows some fundamental unsatis- 
factory caveats of the gapfuU excitation which are not in 
BdG, we believe that the present study is of still some 
use, because the overall features of the resulting vortex 
picture are expected to be independent of the approx- 
imations. These physical quantities obtained here may 
be directly observable once the experiment succeeds in 
creating a vortex. _ 

In his series of papers Fettem^ investigates the vor- 
tex structure in an imperfect Bose gas within the Bo- 
goliubov approximation, giving an approximate analytic 
treatment for an infinite system: Although he succeeded 
in deriving asymptotic behaviors of several physical quan- 
tities, such as the non-condensate fraction at the core 
using the approximately obtained eigenfunctions jOn the 
basis of the work by PitaevskiiQ^ and lordanskiiEP who 
solve the Bogoliubov equation for the particular angular 
momentum states: qg = —1 and qg = 1 respectively {qe 
is defined shortly). Hence the problem with the rotating 
BEC in the Bogoliubov theory is not solved completely 
even for an infinite system. Thus no one has analyzed it 
for a finite trapped system. 

In next Section a brief description of the Hartree-Fock 
Bogoliubov theory (HFB) is given. In this process vari- 
ous approximate treatments mentioned above are intro- 
duced. We explain how to solve a set of the HFB equa- 
tions numerically in a self-consistent mannp;, following 
the method devised by Gygi and Schliiter.113' Here we 
consider a cylindrically symmetric case with a rigid wall 
and choose parameters appropriate to present BEC sys- 
tems for ^^Na and ^^Rb. The results for stationary non- 
uniform case are presented in §3. The vortex structure 
is discussed in details within HFB which turns out to be 
only stable theory among the above theories except for 
CP when the system rotates around the symmetry axis 
in §4. The final section is devoted to discussions and 
summary. 

§2. Formulation and Numerical Procedure 

2.1 Various Approximations 

We start with the following Hamiltonian in which Bose 
particles interact with a two-body potential Vlnt^ 



+ Kxt(r)-Ai^*(r) 



-^jAv^jAv^i>\v^)i>\v2)V,^,{v^-v^)i>{v2)i>{A^) 



where the chemical potential ji is introduced to fix the 
particle number and T4xt(r) is the confining potential. 
In order to describe the Bose condensation, we assume 
that the field operator ^ is decomposed into 



ii = il}(r) + 0(r) 
where the ground state average is given by 
(*(r)) = 0(r). 



(2) 



(3) 



A c-number (l)(r) corresponds to the condensate wave 
function and ip{r) is a q- number describing the non- 
condensate. The two-body interaction Vint(ri ~ r2) is 
assumed to be gS{ri — r2) with g being a positive (re- 
pulsive) constant proportional to the s-wave scattering 
length a, namely g — Anfi^a/m (m the particle mass). 
Substituting the above decomposition (||) in (1), we ob- 
tain 

r(r)/i(r)0(r)-Hi5l0l' 



H 



dr 

V'^(r){;i(r)+.g|0(r)|2}0(r) + 
V>(r){/i(r) + 2.g|0(r)p}V'(r) 
|^t(r)v}t(r)</,(r)</,(r)+h.c. 
g'0^(r)-0(r)'0(r)(/)* (r) + h.c. 



-|^t(r)v}t(r)v3(r)^(r) 



h.c. 



with 



h{r) 



2y72 



rv 

2m 



+ Vcxt(r) - ^J■ 



(4) 



(5) 



H = / dr 



one-body Hamiltonian. Let us introduce the variational 
parameters: the non-condensate density p(r) = {tp^ {r)tp (r)) 
and the anomalous correlation A(r) = (7/)(r)-0(r)) and 
approximate as 

rp'^r/jij! ^2-il!p + il;''A (6) 

'tl)^4!^-iij^ = ^^^^A + 'il>'ipA* +'i'ip^p. (7) 

Then, H is rewritten as 

]r(r)Mr)0(r) + i9|0|4} 

+V?(r){/i(r) + .g|0(r)p -f 2gp{v)}<j>{r) + h.c. 

+-4;'^{r)gA{r)(f)*{r) + h.c. 

+^t(r) {h{r) + 2.9|^(r)|2 + 2gp{v)} ^(r) 

+ {A(r) + ^^(^r)} ^t(r) + h.c] . (8) 

In order to diagonalize this Hamiltonian, the following 
Bogoliubov transformation is employed, namely, ipir) is 
written in terms of the creation and annihilation oper- 
ators rjq and ?7| and the non-condensate wave functions 
Uq{r) and ^^(r) as 



(9) 
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where q denotes the quantum number. This leads to the 
diagonaUzed form: 



(10) 



The condition that the first order term in V'(r) vanish 
yields 

{/i(r)+g|0(r)|2+2.9p(r)}0(r)+gA(r)0*(r) =0. (11) 

When p(r) and A(r) are made zero, it reduces to the 
Gross-Pitaevskii (GP) equation: 



Mr)<^(r)+g0*(r)0(r)0(r) =0, 



(12) 



which is a non-linear Schrodinger type equation. 

The condition that the Hamiltonian be diagonalized 
gives rise to the following set of eigenvalue equations for 
Ug(r) and Wg(r) with the eigenvalue Eq-. 

[Mr) + 25Mr) + |0(r)|2}] «,(r) 

-g [A(r) + 02(r)] Vq{v) = EqUg{r) (13) 

[Mr) + 2.g{p(r) + |0(r)|2}] ^;,(r) 

-5 [A*(r) + 0*^(r)] ^..(r) = ^EgVq{r). (14) 

The eigenvectors Uq{r) and Vq{r) must satisfy the nor- 
malization condition: 

{u*p{r)uq{r) ~ v;{r)vq{r)} dr = Sp^q. (15) 

The variational parameters p{r) and A(r) are determined 
self-consistently by 

p{r) = (tp^) 

= J2i\uqir)f + \vqirr)f{Eq) + J2 Mr)f 
q 1 

A(r) = (^^) 

where f{E) is the Bose distribution function. At abso- 
lute zero temperature, which we consider from now on, 
these become 



q 



(16) 
(17) 



Equations (|l|), (13), (14), (16) and (17) constitute a 
complete set of the self-consistent equations for the Hatree- 
Fock Bogoliubov theory (HFB). If p(r) and A(r) are 
made to zero, we recover the Bogoliubov theory (BA). 
The Popov approximation (PA) is the case of A(r) =0 
in HFB. In the foUowings, we perform these four cases, 
namely, GP, BA, PA and HFB in the equal footing to 
comparatively study similarities and differences. 

The expectation value of the particle number density 
is given as 



that is, the total density consists of the condensate part 
4>{r) and the non-condensate part p{r). The particle cur- 
rent density is calculated, 

j(r) = -^{0*(r)V0(r) - 0(r)Vr (r)} 



2mi 

h 

2mi 



{rP\r) ■ VV'(r) - VV3^(r) • Vi(r)). (19) 



2.2 Numerical Procedure 

For later convenience, we introduce the following non- 
dimensional quantities: In terms of the mass m and the 
average particle number density ria, the various densities 
and the length are scaled by ria and = 
respectively. is the coherence length of the condensate 



given by solving the above GP equation (12). The en- 
ergy is scaled by Uag- We define the quantities: r' = ^r, 
0'(r') ^ ^0(r), p'(r') ^ ^p{v), A'(r') ^ ;i-A(°r), 
u'(r') ^ ^u(r), v'ir') ^ ^v{r), E'q ee ^^Eq. The 
system is now characterized by ria^a ^'^d the system vol- 
ume normalized by . From now on, we suppress primes 
/ in these newly defined quantities. Note that by increas- 
ing the average density Ua the effective interaction gets 
stronger, namely the tuning of the interaction by con- 
trolling an external parameter Ua is an interesting and 
important aspect of the BEG system. 

We now consider a cylindrically symmetric system which 
is characterized by the radius R and the height L. We im- 
pose the boundary conditions that in terms of the cylin- 
drical coordinate: r = (r, 9, z) the all wave functions 
vanish at the wall r ~ R and the periodic boundary 
condition along the z-axis. When a vortex line passes 
through the center of the cylinder the condensate wave 
function (j)(r) is expresses as 



(r, 6, z) — (f>{r)e 



iw9 



(20) 



where 4'{r) is a real function and w is the winding num- 
ber, w — Q corresponds to non-vortex case and w ^ 1 to 
the vortex case. The w > 2 case is not considered here 
because this state is energetically unstable. The non- 
condensate density p is a real function, depending only 
on r that is, p{r,d,z) = p{r). It is seen from eq. ( |l l| ) 
that the anomalous correlation A has the phase with 26*, 
thus A(r, 0, z) — A(r)e^™^. It is also seen from eqs. (13) 
and (14) that the phases of Wg(r) and Wg(r) are written 
as 



(21) 
(22) 



Mq(r) =iiq(r)p>^e'(9<'+"')s 
wq(r) = vq(r)e'«^"e'(*^'")^ 

The set of the quantum numbers q in (|^) is described 
by ((7r, ge, qz) where qr = 1, 2, 3, • • •, qe = 0, ±1, ±2, • • •, 
=0,±27r/L,±4Vi,---. 

The functions Uq{r) and Vq{r) are expanded in terms 

of 

x/2 



Ji 



ImI+1 



Ji, 



^ RJ 



(23) 



(h(r)) = |</)(r)p+p(r). 



(18) 
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as 



i 



(24) 
(25) 

(i) 



where (r) is the Bessel function of /i-th order and a/^ 
denotes i-th zero of J^. 
Then, eq. (^ becomes 

■ d^^ _^ 1 dcj) 
dr^ 7' dr 

+ {K=xt - + 0' + 2p + A} = 0. (26) 

The eigenvalue problem of eqs. (13) and (14) reduces to 
diagonalizing the matrix: 



^Aij{qg +w,qz) 



-Bi^j{qe,w) 



j(2) 



= Ec 



"q 
"q 

V : y 

(27) 




+ (Kxt + 2(02+p))^(;)^OVdr(28) 



B.,, (g,, 7«) = / (A + 02) ^W^^¥,^^;_^r dr (29) 
Jo 

where q = {qr , qe , qz) . The normalization condition (|lq 
of the eigenvectors is rewritten in terms of and tfq 



as 



-^^a E {^P - < } - ^^P.q- (30) 
From eqs. (16) and (17) p{r) and A(r) are determined 

by 



(31) 



A(-) = - E ^S^^qV^lV^.V'^f-^- (32) 



q,i,i 



The iterative calculations of eqs. (26), (p7|), (31), and 
(32) yield a convergent self-consistent solution. It is not 
self-evident a priori that eq. ( p7|) gives real eigenvalues 
because the Hamiltonian matrix in eq. (|2^) is not sym- 
metric. Note that since qe is the angular momentum 
and qz = 2Tm/L (n integer) is the wave number along 
the 0-axis, both being good quantum numbers, the eigen- 
value equation (^ is decomposed into each qg and q^- 
Each block-diagonal eigenvalue equation gives rise to the 
quantum number qr along the radial direction r. 



2.3 Calculated System 

As mentioned, our system is characterized by the two 
parameters: the reduced density 7^aCa ^^.d the system 
size R and L (in unites of ^a). The interaction strength 
characterized by the scattering length a is absorbed in 
the reduced density. We set Ua^^ = l-O' ^ = 1^, and 
L — ?>Q throughout this paper. We are treating 2.12x lO'' 
atoms. Since the scattering length of Rb is known as 
a=5.8nm, the energy scale Uag ~ 10~^°J. These num- 
bers are compared with the experiment by Anderson et 
alH The radius of the condensate ^ lO/xm, the trap- 
ping potential energy ~ 10~2^J and the total number of 
atoms ~10^. Similar set of the numbers are obtained for 
the Na experiment fl) 

We have employed the two methods for numerical cal- 
culation: the momentum cutoff method and the energy 
cutoff method to check our numerical accuracy. In the 
former, the computation is limited to certain finite quan- 
tum numbers for q^, ?6i, and qz. We have examined sev- 
eral cases where -50<g6i, g2<+50 and -75<ge, q^<-l-75, 
keeping the maximum number of qr ~ 50 fixed. In 
view of the tendency of the obtained solutions, the last 
case is satisfactory, thus the total number of the treated 
eigenfunctions ^ 500000 because the size dependence on 
qe and qz almost ceases stopping. In the energy cut- 
off method where the calculation terminates when the 
obtained eigenvalues exceeds a certain value, typically, 
100, and 200 in the units of riag. The last case which 
treats ~ 500000 eigenfunctions seems best in view of 
the expected behavior of the solutions. Note that the 
anomalous correlation A is most dependent on this cut- 
off condition while p is relatively insensitive. Comparing 
the two methods, the energy cutoff method is better than 
the momentum cutoff method, thus we show the results 
employed this method below. 

§3. Stationary non-uniform case 

In order to examine the accuracy of the numerical com- 
putations and to know the properties of BEC in the sta- 
tionary state, we first consider non-vortex Bose system 
confined in a cylindrical vessel. The corresponding infi- 
nite system is analyzed by, for example, Fettem^ within 
the Bogoliubov approximation (BA). 

The condensate (/)(r) is shown in Fig. 1 where GP and 
BA give the same result for this quantity, since BA ne- 
glects the non-condensate and the anomalous correlation. 
It is seen that PA and HFB give nearly same results as 
that in GP (or BA) because the absolute values of p and 
A are very small for our parameter selected. Thus this 
is not the general conclusion. The condensate fraction 
could decrease as g or Ua increases. The overall behav- 
iors of these results are quite understandable; The con- 
densate changes only near the wall, whose characteristic 
length is the coherent length introduced before. The 
expected flatness in the bulk region around the center 
indicates the reliability of our numerical calculations. 

The spatial profiles of the non-condensate p(r) for var- 
ious approximations are displayed in Fig. 2. p(r) in BA 
is calculated by using eq. (16) where u{r) and v{r) are 
a solution in BA and is shown for reference. While BA 
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Fig. 1. Spatial variation of the condensate <f>{r) as a function of 
r for HFB (bold line), PA (dot dashed line), and GP and BA 
(dashed line). 



which is contrasted with other quantities such as (p^r) or 
p{r). It should be noticed also that as mentioned above 
the expected flat behavior far from the wall is reproduced 
only when the enough number of the eigenfunctions are 
taken account in the numerical computation, otherwise 
this particular quantity A(r) often fails to exhibit the 
expected flat behavior. 








-0.3 



and PA give a similar variation, the magnitude of p(r) 
in HFB is almost halved. This difference may be related 
to A in HFB. It is noted that p{0) ^ 0.014 in BA and 
PA compared with the analytic expression at T—Q 

P - , . 0.012^ (33) 



no 



by FetterEf for an infinite system (ng the average number 
density and ^ — Ti / ^2mnog) whose number p = 0.012 in 
the present case (^■^no=l-0). The small difference comes 
from the system size (finite vs infinite). In fact in our 
calculations the limiting value p = 0.012 are slowly re- 
covered as the system size increases. The major spatial 
variation only occurs at the wall whose length is an or- 
der of . Note from Fig. 2 that p vanishes quadratically 
instead of linearly in (/)(r). 



0.005 




Fig. 2. Spatial variation of the non-condensate p(r) as a function 
of r for HFB (bold line), PA (dot dashed line), and BA (dashed 
line) . 



The anomalous correlation A(r) is shown in Fig. 3 
where BA and HFB are compared. The amplitude A(r) 
in HFB at the center are rather large and negative, which 
affects mainly on p{r), namely, the large difference of p(r) 
between BA and PA, and HFB comes from the absence or 
presence of A(r). It is noted also that |A(r)| in the bulk 
region sensitively depends on the calculated system size 




5 Radius 1 



Fig. 3. Spatial variation of the anomalous correlation A(r) as a 
function of r for HFB (bold line) and BA (dashed line). 



In Figs. 4(a), (b), and (c) where the eigenvalues i?(q) 
are plotted as functions of the quantum numbers of qe 
and q^, the excitation spectra are shown for BA and 
HFB. The spectrum in PA which is not shown here is 
almost identical to that in BA. Indeed, the excitation 
spectrumJii PA and BA is gapless while HFB is gapful as 
expected.113) In BA the so-called Bogoliubov spectruim' 
E^{k) is known to be given by 



i?^(k) 



2eknog 



(34) 



where Ck 
to 



In the long wavelength limit it reduces 



£;^(k) - hck for k 







(35) 



where c = \f^^ ■ The low-lying excitations are known 
to be exhausted by phonons (c is the sound velocity). Ej^ 
In the short wavelength limit it becomes 



i?^(k) 



2m 



for k 



(36) 



where the excitations are individual one-particle excita- 
tion. These expected behaviors in BA are well repro- 
duced by the present calculations; As seen from Figs. 
4(a) and 4(b), the dispersion relations in BA and PA are 
linear in g and quadratic in (7 — > cx) where k^ = ^Iz- 
In fact, the theoretical lines of eqs. (35) and (36) drawn 
in Figs. 4(a) and 4(b) show good fits to the numerical 
results without any fitting parameter, proving the reli- 
ability of our calculations. On the other hand, as seen 
from Fig. 4(c) HFB is a gapful theory. This dispersion 
relation is quadratic in both qo and qz. The magnitude 
of this gap is an order of Uag. 
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Fig. 4. Excitation spectra as functions of qg (left hand side) and 
qz (right hand side) in BA for the small wave numbers (a) and the 
overall wave numbers at = (b). The analytical expression 
eq. (35) of the Bogoliubov spectrum in the long-wave length limit 
(straight line) is drawn in the right hand figure of 4(a), and eq. 
(36) in the short wave length limit (lower edge parabola) drawn 
in (b), showing that they agree in those limits, (c) Excitation 
spectra as functions of qg (left hand side) and qz (right hand 
side) in HFB. Note that it has a gap whose magnitude is about 
0.9. 



In Fig. 5 we display the local density of states: 
N{E, r) = 5]{|«q(r)p + kq(r)p}'5(£;q - E) (37) 

q 

a combination of the quasi-particle eigenfunctions with 
the lowest energies in BA (the results for the other ap- 
proximations are quite similar). It is seen from this that 
although the weight spreads out to the entire regions, 
the major one concentrates near the wall boundary. 




15^ 



Fig. 5. Local density of states N{E, r) in the lower energy region 
for BA. It is noted that the density of states accumulated at the 
center r = and the wall r = 15. 



§4. Rotating vortex case 

Having established the reliability of our numerical method, 
we now discuss the results of the isolated singly quan- 
tized vortex case whose quantization unit is ^. It is 
expected that a rotating BEC system whose frequency 
exceeds a certain critical value lUc ^ — ^ In -§ estimated 

[—1 mi?-^ { 

as ^ 50rad/s for ^^Rbc3^ sustains a quantized vortex line 
threading along the cylindrically symmetric z-axis. 

It quickly becomes clear after several numerical trials 
that the BA and PA do not fulfill a fundamental re- 
quirement, that is, the quasi-particle eigenvalues E{q) 
in eqs. (13) and (14) or eq. (27) for BA (p = and 
A = 0) and also for PA (A = 0) must be positive be- 
cause the condensate situates at zero energy. The nega- 
tive eigenvalue means an instability of the vortex state. 
Thus these approximations cannot be a consistent the- 
ory for describing the vortex state. As seen from Ta- 
ble 1 where the quasi-particle eigenvalues £'(q) from the 
lowest energies are listed for BA and HFB, the lowest 
eigenvalue with qg = —1, = and = 1 in BA is 
negative. We have checked the negative eigenvalue in 
BA by changing several conditions: The system size {R, 
L), the interaction strength, and the Hamiltonian matrix 
size. The negative eigenvalue belonging to the lowest en- 
ergy in BA always exists and is not an artifact of our 
numerical computations. As for PA, the lowest eigen- 
value with the same quantum number mentioned above 
becomes also a negative number, which always appears at 
every steps of the iteration processes for self-consistency. 
We never complete a self-consistency and cannot obtain 
a self-consistent solution. We conclude that PA cannot 
sustain a stable vortex solution. Thus, we are left with 
the non self-consistent GP and the full self-consistent 
HFB which are discussed in full details in the following. 

The spatial variations of the condensate 0(r) in HFB 
and GP are shown in Fig. 6. It is seen from this that both 
are almost identical, but the condensate is pushed out in 
HFB by the presence of the non-condensate, resulting in 
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Fig. 6. Spatial variation of the condensate (/!>(r) as a function of r 
for HFB (bold line) and CP (dashed line). It is seen that at the 
vortex core r = it vanishes linearly, recovering its bulk value 
within the characteristic length ^a- 



a slightly larger core radius. This could be further am- 
plified when the interaction g becomes stronger or the 
atomic density ria becomes high since the non-condensate 
fraction increases. It is also noticed that understandably 
4>{r) is almost symmetric at the middle r = 7.5 because 
the characteristic length scale is in this system, which 
governs the spatial variation at the core and wall. At the 
vortex core (j)(0) — and linearly rises to recover its bulk 
value shown in Fig. 1. It will be interesting to check 
if there is a similar Kramer-Pesch effecttiP seen in su- 
perconductors where the core radius shrinks as tempera- 
ture decreases. If indeed exists, the core radius increases 
further as temperature rises. This temperature effect 
belongs to a future problem. This should be checked 
experimentally in the present BEC systems because the 
expected core radius is far larger than that in superfluid 
^He. 




Fig. 7. Spatial variation of the non-condensate p(r) as a func- 
tion of r for HFB. It is seen that p(0) is non-vanishing, coming 
from non-vanishing contribution with qg = 0. The inset shows 
the contributions from various components characterized by the 
quantum numbers: qg = 0, ±1, and ±2. Note the different scales 
of the vertical axes. 



The non-condensate p{r) is displayed in Fig. 7. Ac- 



cording to Fetter,® a universal relation p(0)/p(oo) ~1.4 
at T = is derived for BA, independent of the inter- 
action strength g, that is, the amplitude of p(0) at the 
core must exceed that in the outside region. This pre- 
diction is not supported by the present calculation in 
HFB. On the contrary, our result does show a suppres- 
sion of p{r) around the core region whose characteristic 
length ~ ^a- We do not consider the origin of this dis- 
crepancy further because BA is not a stable theory for 
describing the vortex as mentioned before. In the HFB 
result, which is a stable solution, p{r) recovers the bulk 
value^O.006 (see Fig. 2) far from the core where p(0) re- 
duces to almost zero. The characteristic recovery length 
is evidently longer than that in the condensate as seen 
from Fig. 6. This is partly because the behavior in p(r) 
near the core is quadratic in r while that in 0(r) is linear. 
The main contribution to the non-vanishing p{0) comes 
from the component with qg = 1 as seen from the inset of 
Fig. 7 where other dominant components near the vortex 
core are also depicted. This also explains the quadratic 
behavior in r at the core mathematically. 
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Fig. 8. Spatial variation of the amomalous correlation A(r) as 
a function of r for HFB. The inset shows the contributions 
from various components characterized by the quantum num- 
bers: qe = 0, ±1, and ±2. Note the different scales of the vertial 
axes. 



The anomalous correlation A(r) is shown in Fig. 8 for 
HFB, which is no prior prediction and evaluated here 
for the first time. The sign is negative same as in the 
non-vortex case of previous section and A{r) vanishes 
quadratically at the core and also at the wall, recovering 
its value (~ —0.3) in the bulk (see Fig. 3 for compar- 
ison). As shown as the inset where the contributions 
with the smaller qg, even near the vortex core there are 
no distinctive and/or dominant contributions for A(r). 
As mentioned in §3 this quantity strongly depends on 
the energy cutoff chosen. If the cutoff energy increases, 
|A(r)| grows. This feature is absent in the other quanti- 
ties discussed here. 

In Fig. 9 we exhibit the excitation spectra as functions 
of qg and qz where the distinctive excitations at qg — —1 
are seen, which are isolated from the rest of the con- 
tinuum seen before (Fig. 4(c)). This particular isolated 
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Fig. 9. Excitation spectra as functions of qg (left hand side) and 
Qz (right hand side) in HFB. The distinctive modes at qg = — 1 
corresponds to the Kelvin rwave whose dispersion relation eq. 
(38) obtained by Pitaevskifl) in the long wave length is drawn 
in the right hand side by adjusting the vertical axis for com- 
parison. The bulk of the continuum excitation spectra coincides 
with those in Fig. 4(c). 



excitation is known as the Kelvin wave which is present 
in a vortex of classical liquiiLand corresponds to a heli- 
cal mode of the vortex line.ll3 According to PitaevkiiEP 
who found it in BA, at the long wave length limit the 
dispersion relation is given by 

E{qg = -l,k,) = —^\n— (fc,ea<l). (38) 

Zm Kz^a 

Since this expression is valid only for an infinite system, 
it is hard to judge whether or not our numerical result for 
a finite system agrees with this. Apparently, while the cq. 
(38) is gapless, the present result has a gap. Apart from 
the gap, our result for the dispersion relation does not 
contradict this behavior (see the line of eq. (38) drawn 
Fig. 9). The bulk of the gapful continuum in Fig. 9 is just 
the same as in Fig. 4(c) as expected. 

The local density of states given by eq. (37) with the 
lower energy side are shown in Fig. 10 where the states 
distinctively localized near the core correspond to the 
angular momentum qg = —1. The detailed analyses 
of eqs. (13) and (14) for qej-.= —1 and qg = 1 in BA 
are performed by PitaevskiiJ) and lordanskiitP respec- 
tively: Uq,,g^{r) ~ rl9«+i|(l 0(r2)) and Vq,^g^{r) ~ 
j-ke-il (1 -l_ (9(7.2^) -^vl^ic];! are easily derived by analyzing 
eqs. (13) and (14) for BA. Note that only Mge=-i g2 (0) and 
^98=1,9^(0) are non- vanishing at the core. These proper- 
ties are also true for the full self-consistent HFB. The 
local density of states in superconductop.is directly ob- 
served by scanning tunneling microscopelllP and analyzed 
theoretically within the similar theoretical frameworkM 
quite successfully. Since as mentioned before, the core 
radius in BEC systems is relatively large, there is good 
chance to directly observe the local excitation spectrum 
possibly by an optical method once the vortex can be 
created. 



Finally the circulating current density jg (r) for the 6- 
component which is expressed as 

Jeir)=jl'\r)+ji'\r) (39) 




Fig. 10. Local density of states N{Eq,r) in the lower energy re- 
gion for HFB. It is noted that the density of states accumulates 
at the center r = and the wall r = 15. The distinctive peak 
structures come from those with the quantum number qg = —I, 
indicating that these modes localize at the vortex core. 
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Fig. 11. Current distribution jg (r) which comes from the non- 
condensate where the contributions from the components with 
qg = — 0,±1, and ±2 are shown. The inset shows the total 
current jg{r). Note that the dominant contribution comes from 
the condensate. The current is normalized by — ■ 




where the total current consists of the condensate com- 
ponent J^Q^ (r) and the non-condensate J'q^ (j). The non- 
condensate contribution jf^r) normalized by ^ • ^ 
is depicted for HFB in Fig. 11 where each with qg = 
0, ±1, ±2 is shown separately, and those with smaller \qg\ 
dominate near the core. As is seen from Fig. 11 jg (r) 
in the immediate vicinity of the core is governed by the 
qg ^ component. The negative (positive) qg^s give rise 
to the positive (negative) contribution to jg(r). The inset 
shows the total current density jg (r) where jg (r) ^ r for 
small r and jg{r) ^ ^ for larger r. The relative weights 
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of jg^^ (r) and jg^-* (r) depends on the interaction strength 
and the average density 

§5. Conclusion and Discussions 

We have investigated various approximate mean-field 

type theories (Gross-Pitacvskii, BogoHubov theory, Popov 
theory and Hartree-Fock Bogoliubov theory) within the 
framework of the Bogoliubov approximation for a dilute 
Bosc gas, on which renewed interest is focused recently 
by the discovery of Bose-Einstein condensation in alkali 
atom gases. The above four type theories are numeri- 
cally solved for parameters appropriate to ongoing exper- 
iments on ^^Na and *^Rb atoms and analyzed on an equal 
footing for the first time. We extract several properties 
of experimental interest in BEC systems, namely, the 
spatial structures of the condensate, non-condensate and 
anomalous correlation both in stationary non-uniform 
case and the vortex case binder rotation. A numerical 
procedure for solving these mean-field equations are pre- 
sented and critically assessed for future use. 

In the case of the stationary non-uniform Bosc gas 
confined in a cylindrically symmetric vessel the above 
four theories yield almost identical results for the spatial 
profile of the condensate. The non self-consistent BA 
and self-consistent PA and HFB give similar profiles for 
the non-condensate, but in the last the magnitude are 
halved. 

In the vortex case these mean-field theories are numer- 
ically examined. It is found that GP and PA do not fulfill 
the fundamental requirement, showing an instability of 
the theories, and thus are inadequate for describing a 
vortex. The full self-consistent solution for HFB is ob- 
tained and analyzed in detail. The spatial structures of 
the vortex core for Bose systems; the condensate, non- 
condensate and anomalous correlation are explicitly de- 
rived for the first time. Some characteristics of the local 
density of states and circulating current are pointed out 
in the hope to be observed in BEC systems in alkali atom 
gases. 
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Table I. The lowest energies for the Bogoliubov approximation 
and the Hartree-Fock Bogoliubov theory (HFB). 



